
function y = difD2(w,P,X,Z,R,C)

PZZ = calculatePHat(w,X,Z,R,C);

y = 0;

for ii=1:R
    for jj=1:C 
        
        y = y + myfunc(ii,jj,P,PZZ,C);        
            
    end
end

end



function v = myfunc(ii,jj,P,PZ,C)

p = P(ii,jj);
pHat = PZ(ii,jj);


if jj ~= C
    v = (ii -jj)*(p - pHat);
else
    v = C^2 * (ii-jj)*(p-pHat);
end

end

